1 Introducción

Los datos que se van a analizar en este proyecto han sido obtenidos desde Kaggle. Contienen precios de casas que fueron vendidas desde mayo de 2014 hasta mayo de 2015 en King County que es un condado ubicado en el estado estadounidense de Washington.

2 Objetivo del estudio

Lo que queremos hacer con estos datos es predecir el precio de las casas dependiendo de los datos recogidos.

3 Datos

3.1 Categorización del precio

En nuestro estudio inicial, la variable “Precio” es continua, por lo que vamos a categorizarla. Se van a realizar dos tipos de cetegorizaciones:

  • Categorización_1: se ha categorizado en dos grupos, B1, casas baratas (casas con un precio inferior a 500.000) y casas caras, *C1 (precio> 500.000)
  • Categorización2: se ha categorizado en tres grupos, B2, casas baratas (casas con un precio < 330.000), casas precio medio, M2 (330.000<precio< 650.000) y casas caras, *C2 (precio> 650.000

Para decidir las categorizaciones se ha optado por los cuantiles.

#Categorizamos la variable respuesta price:
quantile(datos$price, prob=seq(0, 1, length = 5))
##      0%     25%     50%     75%    100% 
##   75000  321950  450000  645000 7700000
datos$price_categ1 <- cut(datos$price, breaks = c(0, 500000, 100000000), labels = c("B1", "C1"))
table(datos$price_categ1)
## 
##    B1    C1 
## 12560  9053
datos$price_categ2 <- cut(datos$price, breaks = c(0, 330000, 650000, 100000000), labels = c("B2","M2", "C2"))
table(datos$price_categ2)
## 
##    B2    M2    C2 
##  5881 10525  5207

A continuación se muestran cómo están categorizados los datos:

En este mapa se puede visualizar cómo se distribuyen las casas “Caras” y “Baratas”. Se observa que las casas que están cercanas al agua y a cerca de Seattle por la parte Norte, son más caras y hacia el sur más baratas.

center_lon = median(datos$long,na.rm = TRUE)
center_lat = median(datos$lat,na.rm = TRUE)

factpal2 <- colorFactor(c("green","red"), 
                       datos$price_categ1 )

leaflet(datos) %>% addProviderTiles("Esri.NatGeoWorldMap") %>%
  addCircles(lng = ~long, lat = ~lat,
             color = ~factpal2(datos$price_categ1))  %>%
  # controls
  setView(lng=center_lon, lat=center_lat,zoom = 12) %>%

  addLegend("bottomright", pal = factpal2 , values = ~datos$price_categ1,
            title = "Tipos de Casas",
            opacity = 1)

En este otro mapa se puede visualizar cómo se distribuyen las casas según el precio en tres categorías:“Caras”, Medio y “Baratas”.

center_lon = median(datos$long,na.rm = TRUE)
center_lat = median(datos$lat,na.rm = TRUE)

factpal2 <- colorFactor(c("green","red","yellow"), 
                       datos$price_categ2 )

leaflet(datos) %>% addProviderTiles("Esri.NatGeoWorldMap") %>%
  addCircles(lng = ~long, lat = ~lat,
             color = ~factpal2(datos$price_categ2))  %>%
  # controls
  setView(lng=center_lon, lat=center_lat,zoom = 12) %>%

  addLegend("bottomright", pal = factpal2 , values = ~datos$price_categ2,
            title = "Tipos de Casas 3 categorías",
            opacity = 1)

3.2 Train, test y validación

Se va a separar los datos en los 3 conjuntos de datos fundamentales:

  • Conjunto de datos de entrenamiento: en nuestro estudio datos_train, se corresponde con el 70% del total de los datos.
  • Conjunto de datos de validación: en nuestro estudio datos_validacion, se corresponde con el 15% del total de los datos.
  • Conjunto de datos de test: en nuestro estudia datos_test, se corresponde con el 15% del total de los datos.
num_total=nrow(datos)
set.seed(122556) #reproductividad

# 70% para train
indices_train = sample(1:num_total, .7*num_total)
datos_train = datos[indices_train,]

# 15% para test
indices=seq(1:num_total)
indices_test=indices[-indices_train]
indices_test1 = sample(indices_test, .15*num_total)
datos_test = datos[indices_test1,]

# 15% para validacion
indices_validacion=indices[c(-indices_train,-indices_test1)]
datos_validacion=datos[indices_validacion,]

3.3 Análisis exploratorio

Se van a realizar transformaciones de un conjunto de variables, estas transformaciones se aplicarán a cada conjunto de datos, train, test y validación:

  • Se realiza una transformación logarítmica sobre las variables sqft_living (pies cuadrados de la casa), sqft_lot (pies cuadrados del jardín) y sqft_above (pies cuadrados poe encima del suelo), hay que aclarar que esta última variable esla diferencia entre sqft_living y sqft_basement por lo que va a estar altamentente correlad con sqft_living.

  • Se categorizan las variables:

    • Bathroom, esta varible puede tomar valores decimales de 0.25 en 0.25. El número de baños se contabiliza por las piezas y cada baño completo tiene 4 piezas, por lo que con la nueva agrupación toma valores de 1 a 8 baños.

    • Sqft_basement, se categoriza como 0 las casas que no tienen sótano y 1 las casas que sí tienen sótano.

    • Grade, se va a categorizar del siguiente modo con valor 0=calidad Baja, 1= calidad media y 2= calidad alta.

    • Year_renovated, se categoriza como 0 = no ha tenido renovación y 1 = sí ha tenido renovación.

  • Se pasan a factor las variables: waterfront, view, condition, grade_categ y zipcode

  • Se eliminan Outliers.

3.3.1 Transformaciones datos Train

datos_train <- datos_train[,-2]

datos_train$id <- as.factor(datos_train$id)

datos_train$bathrooms_group <- cut(datos_train$bathrooms,breaks = c(-1,0.25,1,2,3,4,5,6,7,8),labels=c(0,1,2,3,4,5,6,7,8))
datos_train$bathrooms_group <- as.numeric(as.character(datos_train$bathrooms_group))

datos_train$log_sqft_living <- log10(datos_train$sqft_living)
datos_train$log_lot <- log10(datos_train$sqft_lot)
datos_train$log_above <- log10(datos_train$sqft_above)

datos_train$sqft_basement_cat <- cut(datos_train$sqft_basement,breaks = c(-1,0,6000),labels=c(0,1))

datos_train$waterfront<-as.factor(datos_train$waterfront)

datos_train$view<-as.factor(datos_train$view)

datos_train$condition<-as.factor(datos_train$condition)

datos_train$grade_categ <- cut(datos_train$grade, breaks = c(0,4,9,13), labels = c(0,1,2))

datos_train$yr_renovated_catg <-cut(datos_train$yr_renovated, breaks=c(-0.5,1933, 2015), labels= c("0","1"))

datos_train$zipcode<-as.factor(datos_train$zipcode)

3.3.1.1 Eliminación de Outliers

datos_train$posicion<-c(1:nrow(datos_train))
indices_cero_habitaciones<-datos_train[datos_train$bedrooms==0,]$posicion
datos_train<-datos_train[-indices_cero_habitaciones,]

datos_train$posicion<-c(1:nrow(datos_train))
indices_cero_banos<-datos_train$posicion[datos_train$bathrooms_group==0]
datos_train<-datos_train[-indices_cero_banos,]

datos_train$posicion<-c(1:nrow(datos_train))
indice_hab33 <- datos_train[datos_train$bedrooms==33,]$posicion
datos_train[datos_train$posicion == indice_hab33,]$bedrooms = 3

3.3.2 Transformaciones adicionales:

Se han realizado dos categorizaciones adicionales sobre el conjunto de datos. En este caso para ver cómo clasificaban estas variables se ha usado un arbol de decisión, las variables son: zipcode y para bathrooms_group

  • Zipcode, esta variable es de tipo factor y tenía 70 códigos postales, por lo que se ha decidido aplicar un arbol de decisión para ver cómo clasificaba los códigos postales y así volver a categorizarla según el resultado obtenido.
model_selec_zipcode<-rpart(price_categ1~zipcode,data=datos_train ,parms=list(split="gini"))

print(model_selec_zipcode)
## n= 15119 
## 
## node), split, n, loss, yval, (yprob)
##       * denotes terminal node
## 
## 1) root 15119 6385 B1 (0.5776837 0.4223163)  
##   2) zipcode=98001,98002,98003,98010,98011,98014,98019,98022,98023,98024,98028,98030,98031,98032,98034,98038,98042,98045,98055,98056,98058,98059,98070,98092,98106,98108,98118,98125,98126,98133,98144,98146,98148,98155,98166,98168,98178,98188,98198 8243 1336 B1 (0.8379231 0.1620769) *
##   3) zipcode=98004,98005,98006,98007,98008,98027,98029,98033,98039,98040,98052,98053,98065,98072,98074,98075,98077,98102,98103,98105,98107,98109,98112,98115,98116,98117,98119,98122,98136,98177,98199 6876 1827 C1 (0.2657068 0.7342932) *

Se va a categorizar en dos zona1 y Zona2.

datos_train$zona<-recode(datos_train$zipcode, "98001=1; 98002=1; 98003=1; 98010=1; 98011=1; 98014=1; 98019=1; 98022=1; 98023=1; 98024=1; 98028=1; 98030=1; 98031=1; 98032=1; 98034=1; 98038=1; 98042=1; 98045=1; 98055=1; 98056=1; 98058=1; 98059=1; 98070=1; 98092=1 ;98106=1; 98108=1; 98118=1; 98125=1; 98126=1; 98133=1; 98144=1; 98146=1; 98148=1; 98155=1; 98166=1; 98168=1; 98178=1; 98188=1; 98198=1; 98004=2; 98005=2; 98006=2; 98007=2; 98008=2; 98027=2; 98029=2; 98033=2; 98039=2; 98040=2; 98052=2; 98053=2; 98065=2; 98072=2; 98074=2; 98075=2; 98077=2; 98102=2; 98103=2; 98105=2; 98107=2; 98109=2; 98112=2; 98115=2; 98116=2; 98117=2; 98119=2; 98122=2; 98136=2; 98177=2; 98199=2")

datos_train$zipcode = NULL
  • bathrooms_group, se aplica el mismo método que con zipcode para ver cómo se puede categorizar esta variable. Toma valores de 1 a 8 y queremos reducir el número de níveles.
model_selec_bathrooms<-rpart(price_categ1~bathrooms_group,data=datos_train ,parms=list(split="gini"))

print(model_selec_bathrooms)
## n= 15119 
## 
## node), split, n, loss, yval, (yprob)
##       * denotes terminal node
## 
## 1) root 15119 6385 B1 (0.5776837 0.4223163)  
##   2) bathrooms_group< 2.5 7282 1850 B1 (0.7459489 0.2540511) *
##   3) bathrooms_group>=2.5 7837 3302 C1 (0.4213347 0.5786653) *

En el resultado del modelo se ve que corta en el número de baños <2.5, por lo que se va a categorizar como 0 aquellas casas que tengan de 1 a 2 baños y como 1 las casas que tengan más de 2 baños.

datos_train$bathrooms_group <- cut(datos_train$bathrooms_group, breaks = c(-1,2.5,8),labels=c(0,1))
#0 1ó 2 baños, 1 de 3 baños en adelante
datos_train_limpio <- datos_train[c(3,22:26,8:10,27,16,17,30,28,20,21)]

#Eliminamos sqft_above porque es una combinalción lineal de sqft_living, están altamente correladas........

datos_train_limpio$log_above = NULL


datos_train_numeric <- datos_train_limpio %>% select_if(is.numeric)

3.3.3 Transformaciones datos Test

Realizamos todas las transformaciones y categorizaciones para el conjunto de datos de test.

datos_test <- datos_test[,-2]

datos_test$id <- as.factor(datos_test$id)

datos_test$bathrooms_group <- cut(datos_test$bathrooms,breaks = c(-1,0.25,1,2,3,4,5,6,7,8),labels=c(0,1,2,3,4,5,6,7,8))
datos_test$bathrooms_group <- as.numeric(as.character(datos_test$bathrooms_group))

datos_test$log_sqft_living <- log10(datos_test$sqft_living)
datos_test$log_lot <- log10(datos_test$sqft_lot)
datos_test$log_above <- log10(datos_test$sqft_above)

datos_test$sqft_basement_cat <- cut(datos_test$sqft_basement,breaks = c(-1,0,6000),labels=c(0,1))

datos_test$waterfront<-as.factor(datos_test$waterfront)

datos_test$view<-as.factor(datos_test$view)

datos_test$condition<-as.factor(datos_test$condition)

datos_test$grade_categ <- cut(datos_test$grade, breaks = c(0,4,9,13), labels = c(0,1,2))

datos_test$yr_renovated_catg <-cut(datos_test$yr_renovated, breaks=c(-0.5,1933, 2015), labels= c("0","1"))

datos_test$zipcode<-as.factor(datos_test$zipcode)

#codificar la variable Zipcode

datos_test$zona<-recode(datos_test$zipcode, " 98001=1; 98002=1; 98003=1; 98010=1; 98011=1; 98014=1; 98019=1; 98022=1; 98023=1; 98024=1; 98028=1; 98030=1; 98031=1; 98032=1; 98034=1; 98038=1; 98042=1; 98045=1; 98055=1; 98056=1; 98058=1; 98059=1; 98070=1; 98092=1 ;98106=1; 98108=1; 98118=1; 98125=1; 98126=1; 98133=1; 98144=1; 98146=1; 98148=1; 98155=1; 98166=1; 98168=1; 98178=1; 98188=1; 98198=1; 98004=2; 98005=2; 98006=2; 98007=2; 98008=2; 98027=2; 98029=2; 98033=2; 98039=2; 98040=2; 98052=2; 98053=2; 98065=2; 98072=2; 98074=2; 98075=2; 98077=2; 98102=2; 98103=2; 98105=2; 98107=2; 98109=2; 98112=2; 98115=2; 98116=2; 98117=2; 98119=2; 98122=2; 98136=2; 98177=2; 98199=2")

datos_test$zipcode = NULL

datos_test$bathrooms_group <- cut(datos_test$bathrooms_group, breaks = c(-1,2.5,8),labels=c(0,1))

#0 1ó 2 baños, 1 de 3 baños en adelante
datos_test_limpio <- datos_test[c(3,22:26,8:10,27,16,17,29,28,20,21)]
datos_test_limpio$log_above = NULL
datos_test_numeric <- datos_test_limpio %>% select_if(is.numeric)

3.4 Transformaciones datos Validación

Realizamos todas las transformaciones y categorizaciones para el conjunto de datos de Validación.

datos_validacion <- datos_validacion[,-2]

datos_validacion$id <- as.factor(datos_validacion$id)

datos_validacion$bathrooms_group <- cut(datos_validacion$bathrooms,breaks = c(-1,0.25,1,2,3,4,5,6,7,8),labels=c(0,1,2,3,4,5,6,7,8))
datos_validacion$bathrooms_group <- as.numeric(as.character(datos_validacion$bathrooms_group))

datos_validacion$log_sqft_living <- log10(datos_validacion$sqft_living)
datos_validacion$log_lot <- log10(datos_validacion$sqft_lot)
datos_validacion$log_above <- log10(datos_validacion$sqft_above)

datos_validacion$sqft_basement_cat <- cut(datos_validacion$sqft_basement,breaks = c(-1,0,6000),labels=c(0,1))

datos_validacion$waterfront<-as.factor(datos_validacion$waterfront)

datos_validacion$view<-as.factor(datos_validacion$view)

datos_validacion$condition<-as.factor(datos_validacion$condition)

datos_validacion$grade_categ <- cut(datos_validacion$grade, breaks = c(0,4,9,13), labels = c(0,1,2))

datos_validacion$yr_renovated_catg <-cut(datos_validacion$yr_renovated, breaks=c(-0.5,1933, 2015), labels= c("0","1"))

datos_validacion$zipcode<-as.factor(datos_validacion$zipcode)

#codificar la variable Zipcode

datos_validacion$zona<-recode(datos_validacion$zipcode, " 98001=1; 98002=1; 98003=1; 98010=1; 98011=1; 98014=1; 98019=1; 98022=1; 98023=1; 98024=1; 98028=1; 98030=1; 98031=1; 98032=1; 98034=1; 98038=1; 98042=1; 98045=1; 98055=1; 98056=1; 98058=1; 98059=1; 98070=1; 98092=1 ;98106=1; 98108=1; 98118=1; 98125=1; 98126=1; 98133=1; 98144=1; 98146=1; 98148=1; 98155=1; 98166=1; 98168=1; 98178=1; 98188=1; 98198=1; 98004=2; 98005=2; 98006=2; 98007=2; 98008=2; 98027=2; 98029=2; 98033=2; 98039=2; 98040=2; 98052=2; 98053=2; 98065=2; 98072=2; 98074=2; 98075=2; 98077=2; 98102=2; 98103=2; 98105=2; 98107=2; 98109=2; 98112=2; 98115=2; 98116=2; 98117=2; 98119=2; 98122=2; 98136=2; 98177=2; 98199=2")

datos_validacion$zipcode = NULL

datos_validacion$bathrooms_group <- cut(datos_validacion$bathrooms_group, breaks = c(-1,2.5,8),labels=c(0,1))
datos_validacion_limpio <- datos_validacion[c(3,22:26,8:10,27,16,17,29,28,20,21)]
datos_validacion_limpio$log_above = NULL

datos_validacion_numeric <- datos_validacion_limpio %>% select_if(is.numeric)
#Aquí nos cargamos ya price_categ2

datos_train_limpio$price_categ2= NULL
datos_test_limpio$price_categ2= NULL
datos_validacion_limpio$price_categ2 = NULL

4 Métodos de agrupamiento-No supervisado

A continuación, se implementan dos métodos de agrupamiento no supervisado. Primero, un método jerárquico para averiguar (si es posible) la cantidad adecuada de clústeres y a continuación K-Means.

4.1 Cluster Jerárquico

La agrupación es una técnica para agrupar puntos de datos similares en un grupo y separar las diferentes observaciones en diferentes grupos o grupos. En el Clustering Jerárquico los clusters se crean de manera que tengan un orden predeterminado. Existen En nuestro estudio se va a plicar un método aglomerativo que consiste en que cada observación se asigna a su propio clúster. Luego, se calcula la similitud (o distancia) entre cada uno de los clusters y los dos clusters más similares se fusionan en uno. Finalmente, los pasos 2 y 3 se repiten hasta que solo quede un grupo.

Aplicando este método a nuestros datos queremos observar si siguen algún patrón para poder agrupar las casas.

4.1.1 Train

Escalamos las variables numéricas, es decir cada variable ahora tendrá una media cero y una desviación estándar uno.También se requiere los valores de distancia.La medida predeterminada para la función dist es ‘Euclidiana’.

datos_scale<- as.data.frame(scale(datos_train_numeric))

#Matriz de distancias

matriz_dist=dist(datos_scale)

En este caso particular, estimamos dos clústeres:

modelo_hc= hclust(matriz_dist, method = "average")

plot(modelo_hc)
rect.hclust(modelo_hc, k=2,border="red")

A continuación vemos como se han agrupado los datos marcando que el número de clústeres sea 2. Prácticamente todas las casas están en el grupo 1.

grupos2=cutree(modelo_hc,k=2)
table(datos_train_limpio$price_categ1, grupos2)
##     grupos2
##         1    2
##   B1 8724   10
##   C1 6385    0

4.2 Clúster no Jerárquico- K-Means

El método de K-Means basa su funcionamiento en agrupar los datos de entrada en un total de k conjuntos definidos por un centroide, cuya distancia con los puntos que pertenecen a cada uno de los datos es la menor posible.

4.2.1 Train

Primero se van a realizar los gráficos para ver cómo están diferenciadas las casa por precio. Para poder visualizarlo en dos dimensiones se ha usado la función “prcomp” (PCA)

colores1= c("red","blue")
colores11 = colores1[datos_train_limpio$price_categ1]

plot(prcomp(datos_train_numeric, scale = T)$x[,1:2], type="n",main= "Dos categorías")
text(prcomp(datos_train_numeric, scale = T)$x[,1:2], as.character(datos_train_limpio$price_categ1), col=colores11)

A continuación se va a aplicar el método de agrupamiento de K-means, al igual que en el método anterior se le pasa la matriz de distacias y se van a agrupar los datos en dos conjuntos:

set.seed(1234)
modelkm1 <- kmeans(matriz_dist, centers=2)
#modelkm10 <- kmeans(matriz_dist, centers=10)
#table(modelkm10$cluster)



#modelkm2 <- kmeans(matriz_dist, centers=2)

#table(modelkm2$cluster) #asignación de observación a los cluster
# modelkm1$totss  #Inercia total
# modelkm1$betweenss  #Inercia inter grupos
# modelkm1$withinss   #Inercia intra grupos
# modelkm1$tot.withinss  #inercia total intra grupos
colores1= c("red","blue")
colores11 = colores1[datos_train_limpio$price_categ1]

plot(prcomp(datos_train_numeric, scale = T)$x[,1:2], type="n")
text(prcomp(datos_train_numeric, scale = T)$x[,1:2], as.character(modelkm1$cluster), col=colores11)

4.2.1.1 Cálculo del k-óptimo.

Se va a determinar la cantidad óptima de centroides a utilizar a partir del Método del Codo. Para ello, aplicaremos la función kmeans al conjunto de datos, variando en cada caso el valor de k, y acumulando los valores de WCSS obtenidos:

set.seed(1234)
wcss <- vector()
for(i in 1:20){
  wcss[i] <- sum(kmeans(datos_scale, i)$withinss)
}
## Warning: did not converge in 10 iterations
ggplot() + geom_point(aes(x = 1:20, y = wcss), color = 'blue') + 
  geom_line(aes(x = 1:20, y = wcss), color = 'blue') + 
  ggtitle("Método del Codo") + 
  xlab('Cantidad de Centroides k') + 
  ylab('WCSS')

Como se observa en la gráfica, el K-óptimo que se podría aplicar sería de 10.

A continuación para visualizar si el agrupamiento que se ha llevado a cabo está relacionado con el precio de las casas(“Caras”, Baratas), se representa en el mapa que se muestra a continuación.

clusterkmeans=as.data.frame(modelkm1$cluster)

clusterkmeans$indice= as.integer(rownames(clusterkmeans))

colnames(clusterkmeans)[1]= "categoria_price_km"
clustering1= clusterkmeans[order(clusterkmeans$indice),]

Se visualiza como se reparten las casas según el agrupamiento que se ha realizado con K-Means:

center_lon = median(datos_train_limpio$long,na.rm = TRUE)
center_lat = median(datos_train_limpio$lat,na.rm = TRUE)

factpal1 <- colorFactor(c("green","red"),
                       clustering1$categoria_price_km )

leaflet(datos_train_limpio) %>% addProviderTiles("Esri.NatGeoWorldMap") %>%
  addCircles(lng = ~long, lat = ~lat,
             color = ~factpal1(clustering1$categoria_price_km))  %>%
  # controls
  setView(lng=center_lon, lat=center_lat,zoom = 12) %>%

  addLegend("bottomright", pal = factpal1 , values = ~clustering1$categoria_price_km,
            title = "Tipos de casas",
            opacity = 1)

Cómo se observa en el mapa, y si lo comparamos con el de los datos iniciales, dista bastante. Por lo que deducimos que la agrupación que está relalizando K-Means no es muy buena.

4.3 K-Medoids

K-medoids es un método de clustering muy similar a K-means en cuanto a que ambos agrupan las observaciones en K clusters, donde K es un valor preestablecido. La diferencia es que, en K-medoids, cada cluster está representado por una observación presente en el cluster (medoid),en nuestro estudio será una observación de una casa, mientras que en K-means cada cluster está representado por su centroide, que se corresponde con el promedio de todas las observaciones del cluster pero con ninguna en particular.

Este algoritmo es menos sensible al ruido y los valores atípicos, en comparación con k-means, porque usa medoides como centros de clúster en lugar de centroides. (utilizados en k-means).

# datos_train_limpio[,1:3]
# dist_eu <- as.matrix(dist(datos_train_limpio[,c(1,3)]))
# gower_daisy <- as.matrix(daisy(datos_train_limpio[,c(1,3)], metric = 'gower'))
# gower_daisy[1:10, 1:10]
datoskmedoids = datos_train_limpio[,-15]
medoids_model = pam(x = datoskmedoids, k = 2, keep.diss = TRUE, keep.data = TRUE)
medoids_model$medoids
##      bedrooms bathrooms_group log_sqft_living  log_lot sqft_basement_cat
## 9606        3               1        3.152288 3.870404                 1
## 8735        4               2        3.406540 3.936111                 2
##      waterfront view condition grade_categ     lat     long zona
## 9606          1    1         3           2 47.4949 -122.239    1
## 8735          1    1         3           2 47.6477 -122.197    2
##      yr_renovated_catg price_categ1
## 9606                 1            1
## 8735                 1            2
# plot(medoids_model, main = "silhouette plot")
grupos<-data.frame(datoskmedoids)
grupos<-cbind(grupos,data.frame(medoids_model$clustering))
grupos$medoids_model.clustering<-as.factor(grupos$medoids_model.clustering)
ggpairs(grupos,columns=c(1,2,3,15,12),mapping=aes(color=medoids_model.clustering))
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

table(grupos$medoids_model.clustering)
## 
##    1    2 
## 8831 6288
#medoids_model$clustering
clustering= sort(medoids_model$clustering)

clustering=as.data.frame(medoids_model$clustering)

clustering$indice= as.integer(rownames(clustering))

colnames(clustering)[1]= "categoria_price"
clustering2= clustering[order(clustering$indice),]
center_lon = median(datoskmedoids$long,na.rm = TRUE)
center_lat = median(datoskmedoids$lat,na.rm = TRUE)

factpal2 <- colorFactor(c("green","red"), 
                       clustering2$categoria_price )

leaflet(datoskmedoids) %>% addProviderTiles("Esri.NatGeoWorldMap") %>%
  addCircles(lng = ~long, lat = ~lat,
             color = ~factpal2(clustering2$categoria_price))  %>%
  # controls
  setView(lng=center_lon, lat=center_lat,zoom = 12) %>%

  addLegend("bottomright", pal = factpal2 , values = ~clustering2$categoria_price,
            title = "Tipos de casas",
            opacity = 1)

5 Técnicas de reducción de la dimensionalidad:

5.1 PCA

Es un método que permite simplificar la complejidad de espacios muestrales con muchas dimensiones a la vez que conserva su información.

pca<-prcomp(datos_train_numeric,scale=T)
plot(prcomp(datos_train_numeric,scale=T))

#pca$center para saber la media de cada variable
#pca$scale para saber la desv.típica de cada variable
summary(prcomp(datos_train_numeric, scale=T))
## Importance of components:
##                          PC1    PC2    PC3    PC4     PC5
## Standard deviation     1.414 1.0885 0.9290 0.7850 0.57947
## Proportion of Variance 0.400 0.2369 0.1726 0.1232 0.06716
## Cumulative Proportion  0.400 0.6370 0.8096 0.9328 1.00000
biplot(x = pca, scale = 0, cex = 0.6, col = c("blue4", "brown3"))

6 Aprendizaje supervisado

6.1 GLM-Regresión Logística.

Con este modelo se va a estudiar si existe relación entre el hecho de que una casa sea “cara” ó “barata” dependiendo de las características de las casas.Se va a generar un modelo en el que a apartir de las variables,…..prediga la probabilidad de que una casa sea barata o cara.

6.1.1 Train

datos_train_rl <- datos_train_limpio[,-c(8,9)]
datos_train_rl$price_categ1<- recode(datos_train_rl$price_categ1, "'B1'=0; 'C1'=1")


#model_glm = glm(price ~., family = binomial, data =datos_train)
model_glm = glm(price_categ1 ~., family = binomial, data =datos_train_rl )
summary(model_glm)
## 
## Call:
## glm(formula = price_categ1 ~ ., family = binomial, data = datos_train_rl)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -3.3536  -0.3998  -0.0883   0.3481   3.6096  
## 
## Coefficients:
##                      Estimate Std. Error z value Pr(>|z|)    
## (Intercept)        -759.91116   32.36457 -23.480  < 2e-16 ***
## bedrooms             -0.23062    0.03834  -6.015 1.80e-09 ***
## bathrooms_group1      0.14717    0.06828   2.155  0.03113 *  
## log_sqft_living      13.13156    0.32565  40.324  < 2e-16 ***
## log_lot               0.66619    0.08223   8.101 5.45e-16 ***
## sqft_basement_cat1   -0.45846    0.05867  -7.814 5.52e-15 ***
## waterfront1           1.35682    0.62801   2.161  0.03073 *  
## view1                 1.13634    0.22417   5.069 4.00e-07 ***
## view2                 1.21640    0.14137   8.605  < 2e-16 ***
## view3                 1.84563    0.21019   8.781  < 2e-16 ***
## view4                 2.49482    0.49846   5.005 5.58e-07 ***
## lat                   5.32574    0.24476  21.759  < 2e-16 ***
## long                 -3.75931    0.24366 -15.429  < 2e-16 ***
## zona2                 3.33880    0.07036  47.451  < 2e-16 ***
## yr_renovated_catg1    0.37667    0.13792   2.731  0.00631 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 20592.9  on 15118  degrees of freedom
## Residual deviance:  8975.8  on 15104  degrees of freedom
## AIC: 9005.8
## 
## Number of Fisher Scoring iterations: 6
predicciones <- ifelse(test = model_glm$fitted.values > 0.5, yes = 1, no = 0)
tabla_reglog <- table(model_glm$model$price_categ1, predicciones,
                          dnn = c("observaciones", "predicciones"))

#caras

pred.means2=tabla_reglog[2,2]/(tabla_reglog[2,2]+tabla_reglog[1,2])
rec.means2=tabla_reglog[2,2]/(tabla_reglog[2,2]+tabla_reglog[2,1])
F_caras_reglog=(2*pred.means2*rec.means2)/(pred.means2+rec.means2)
F_caras_reglog
## [1] 0.8477249
tabla_reglog
##              predicciones
## observaciones    0    1
##             0 7813  921
##             1 1010 5375
#baratas
pred.means2=tabla_reglog[1,1]/(tabla_reglog[1,1]+tabla_reglog[2,1])
rec.means2=tabla_reglog[1,1]/(tabla_reglog[1,1]+tabla_reglog[1,2])
F_baratas_reglog=(2*pred.means2*rec.means2)/(pred.means2+rec.means2)
F_baratas_reglog
## [1] 0.8900154
tabla_reglog
##              predicciones
## observaciones    0    1
##             0 7813  921
##             1 1010 5375
#F-Medida

F_reglog_train= (F_caras_reglog+F_baratas_reglog)/2
F_reglog_train
## [1] 0.8688702
accuracy_reglog_train = (tabla_reglog[1,1]+tabla_reglog[2,2]) / (tabla_reglog[1,1]+tabla_reglog[1,2]+tabla_reglog[2,1]+tabla_reglog[2,2])
accuracy_reglog_train
## [1] 0.8722799

Explicamos con palabras porqué este modelo no es bueno. En vez de seguir con test y validación se descartan.

#ggplot(subset(datos_train_rl, select = c(bedrooms,bathrooms_group, #log_sqft_living, log_lot, log_above, lat, long, price_categ1)), aes(x = #bathrooms_group,
#y = price_categ1)) + geom_point() +  # geom_smooth(method = 'glm', method.args = list(family =
# 'binomial'), se=FALSE) +
#geom_line(data = data.frame(bathrooms_group = data.frame(bathrooms_group = #seq(1,8,1)), price_categ1 = predict(model_glm, data.frame(bathrooms_group = #seq(1,8,1)), type = "response")), colour = "navy")

6.2 KNN

Ver con cuantos vecinos:

suppressWarnings(suppressMessages(library(kknn)))
modelo <- train.kknn(price_categ1 ~ ., data = datos_train_limpio, kmax = 20)
modelo
## 
## Call:
## train.kknn(formula = price_categ1 ~ ., data = datos_train_limpio,     kmax = 20)
## 
## Type of response variable: nominal
## Minimal misclassification: 0.1148886
## Best kernel: optimal
## Best k: 15

6.2.1 Train

#k=1
#train
constante = 1
knn.train1=knn.cv(scale(datos_train_numeric),cl=as.factor(datos_train_limpio$price_categ1),k=constante)
tabla.knn.train1=table(knn.train1,datos_train_limpio$price_categ1)
pred.knn.train1=tabla.knn.train1[2,2]/(tabla.knn.train1[2,2]+tabla.knn.train1[1,2])
rec.knn.train1=tabla.knn.train1[2,2]/(tabla.knn.train1[2,2]+tabla.knn.train1[2,1])
F_medida.knn.train1=(5*pred.knn.train1*rec.knn.train1)/(4*pred.knn.train1+rec.knn.train1)
F_medida.knn.train1
## [1] 0.8367726
tabla.knn.train1
##           
## knn.train1   B1   C1
##         B1 7707 1073
##         C1 1027 5312
#k=2
#train
knn_train2=knn.cv(scale(datos_train_numeric),cl=as.factor(datos_train_limpio$price_categ1),k=2)
tabla_knn_train2=table(knn_train2,datos_train_limpio$price_categ1)
pred_knn_train2=tabla_knn_train2[2,2]/(tabla_knn_train2[2,2]+tabla_knn_train2[1,2])
rec_knn_train2=tabla_knn_train2[2,2]/(tabla_knn_train2[2,2]+tabla_knn_train2[2,1])
F_medida_knn_train2=(5*pred_knn_train2*rec_knn_train2)/(4*pred_knn_train2+rec_knn_train2)
F_medida_knn_train2
## [1] 0.8304363
tabla_knn_train2
##           
## knn_train2   B1   C1
##         B1 7651 1082
##         C1 1083 5303
center_lon = median(datos_train_limpio$long,na.rm = TRUE)
center_lat = median(datos_train_limpio$lat,na.rm = TRUE)

factpal1 <- colorFactor(c("green","red"), 
                       knn_train2)

leaflet(datos_train_limpio) %>% addProviderTiles("Esri.NatGeoWorldMap") %>%
  addCircles(lng = ~long, lat = ~lat,
             color = ~factpal1(knn_train2))  %>%
  # controls
  setView(lng=center_lon, lat=center_lat,zoom = 12) %>%

  addLegend("bottomright", pal = factpal1 , values = ~knn_train2,
            title = "Precio (en miles de $)",
            opacity = 1)
tabla_knn=table(datos_train_limpio$price_categ1,knn_train2)

#caras
pred_means_caras_knn=tabla_knn[2,2]/(tabla_knn[2,2]+tabla_knn[1,2])
rec_means_knn=tabla_knn[2,2]/(tabla_knn[2,2]+tabla_knn[2,1])
F_caras_knn_train=(2*pred_means_caras_knn*rec_means_knn)/(pred_means_caras_knn+rec_means_knn)
F_caras_knn_train
## [1] 0.8304753
tabla_knn
##     knn_train2
##        B1   C1
##   B1 7651 1083
##   C1 1082 5303
#baratas
pred_means_baratas_knn=tabla_knn[1,1]/(tabla_knn[1,1]+tabla_knn[2,1])
rec_means_baratas_knn=tabla_knn[1,1]/(tabla_knn[1,1]+tabla_knn[1,2])
F_baratas_knn_train=(2*pred_means_baratas_knn*rec_means_baratas_knn)/(pred_means_baratas_knn+rec_means_baratas_knn)
F_baratas_knn_train
## [1] 0.876052
tabla_knn
##     knn_train2
##        B1   C1
##   B1 7651 1083
##   C1 1082 5303
#Media de la F-MEDIDA para KNN
F_knn_train= (F_caras_knn_train+F_baratas_knn_train)/2
F_knn_train
## [1] 0.8532636
accuracy_knn_train= (tabla_knn[1,1]+tabla_knn[2,2]) / (tabla_knn[1,1]+tabla_knn[1,2]+tabla_knn[2,1]+tabla_knn[2,2])
accuracy_knn_train
## [1] 0.8568027

6.2.2 Test

knn_test2=knn(scale(datos_train_numeric),scale(datos_test_numeric),cl=datos_train_limpio$price_categ1,k=2)

tabla2_test=table(datos_test_limpio$price_categ1,knn_test2)
#caras
pred.means2=tabla2_test[2,2]/(tabla2_test[2,2]+tabla2_test[1,2])
rec.means2=tabla2_test[2,2]/(tabla2_test[2,2]+tabla2_test[2,1])
F_medida.means2=(2*pred.means2*rec.means2)/(pred.means2+rec.means2)
F_medida.means2
## [1] 0.8109322
tabla2_test
##     knn_test2
##        B1   C1
##   B1 1653  258
##   C1  247 1083
#baratas
pred.means2=tabla2_test[1,1]/(tabla2_test[1,1]+tabla2_test[2,1])
rec.means2=tabla2_test[1,1]/(tabla2_test[1,1]+tabla2_test[1,2])
F_medida.means2=(2*pred.means2*rec.means2)/(pred.means2+rec.means2)
F_medida.means2
## [1] 0.8674888
tabla2_test
##     knn_test2
##        B1   C1
##   B1 1653  258
##   C1  247 1083
accuracy = (tabla2_test[1,1]+tabla2_test[2,2]) / (tabla2_test[1,1]+tabla2_test[1,2]+tabla2_test[2,1]+tabla2_test[2,2])
accuracy
## [1] 0.8441839

6.2.3 Validación

knn.val=knn(scale(datos_train_numeric),scale(datos_validacion_numeric),cl=datos_train_limpio$price_categ1,k=2)

tabla2_val=table(datos_validacion_limpio$price_categ1,knn.val)
#caras
pred.means2=tabla2_val[2,2]/(tabla2_val[2,2]+tabla2_val[1,2])
rec.means2=tabla2_val[2,2]/(tabla2_val[2,2]+tabla2_val[2,1])
F_medida.means2=(2*pred.means2*rec.means2)/(pred.means2+rec.means2)
F_medida.means2
## [1] 0.8139274
tabla2_val
##     knn.val
##        B1   C1
##   B1 1659  247
##   C1  250 1087
#baratas
pred.means2=tabla2_val[1,1]/(tabla2_val[1,1]+tabla2_val[2,1])
rec.means2=tabla2_val[1,1]/(tabla2_val[1,1]+tabla2_val[1,2])
F_medida.means2=(2*pred.means2*rec.means2)/(pred.means2+rec.means2)
F_medida.means2
## [1] 0.8697248
tabla2_val
##     knn.val
##        B1   C1
##   B1 1659  247
##   C1  250 1087
accuracy = (tabla2_val[1,1]+tabla2_val[2,2]) / (tabla2_val[1,1]+tabla2_val[1,2]+tabla2_val[2,1]+tabla2_val[2,2])
accuracy
## [1] 0.8467468

6.3 Decision Trees

6.3.1 Train

#quitamos lat y long:
datos_decision_tree<-datos_train_limpio[,-c(11,12)]

decisiontree_model=rpart(price_categ1~., data=datos_decision_tree,parms=list(split="gini") )
#resumen=summary(decisiontree_model)
print(decisiontree_model)
## n= 15119 
## 
## node), split, n, loss, yval, (yprob)
##       * denotes terminal node
## 
##  1) root 15119 6385 B1 (0.57768371 0.42231629)  
##    2) log_sqft_living< 3.35698 9914 2474 B1 (0.75045390 0.24954610)  
##      4) lat< 47.55595 4723  252 B1 (0.94664408 0.05335592) *
##      5) lat>=47.55595 5191 2222 B1 (0.57195145 0.42804855)  
##       10) lat>=47.69045 2056  353 B1 (0.82830739 0.17169261) *
##       11) lat< 47.69045 3135 1266 C1 (0.40382775 0.59617225)  
##         22) log_sqft_living< 3.173768 1313  467 B1 (0.64432597 0.35567403) *
##         23) log_sqft_living>=3.173768 1822  420 C1 (0.23051592 0.76948408) *
##    3) log_sqft_living>=3.35698 5205 1294 C1 (0.24860711 0.75139289)  
##      6) lat< 47.51495 1422  451 B1 (0.68284107 0.31715893)  
##       12) log_lot< 4.175105 973  183 B1 (0.81192189 0.18807811) *
##       13) log_lot>=4.175105 449  181 C1 (0.40311804 0.59688196) *
##      7) lat>=47.51495 3783  323 C1 (0.08538197 0.91461803) *
#kk=as.data.frame(summary(decisiontree_model))

fancyRpartPlot(decisiontree_model)

#rpart.plot(decisiontree_model)
#prp(decisiontree_model, type = 2, extra = 6, split.font = 3)
tabla_train_arbol=table(obs = datos_decision_tree$price_categ1, pred = predict(decisiontree_model, type = "class"))

tabla_train_arbol
##     pred
## obs    B1   C1
##   B1 7810  924
##   C1 1255 5130
#caras
pred.means2=tabla_train_arbol[2,2]/(tabla_train_arbol[2,2]+tabla_train_arbol[1,2])
rec.means2=tabla_train_arbol[2,2]/(tabla_train_arbol[2,2]+tabla_train_arbol[2,1])
F_caras_dt_train=(2*pred.means2*rec.means2)/(pred.means2+rec.means2)
F_caras_dt_train
## [1] 0.8248251
tabla_train_arbol
##     pred
## obs    B1   C1
##   B1 7810  924
##   C1 1255 5130
#baratas
pred.means2=tabla_train_arbol[1,1]/(tabla_train_arbol[1,1]+tabla_train_arbol[2,1])
rec.means2=tabla_train_arbol[1,1]/(tabla_train_arbol[1,1]+tabla_train_arbol[1,2])
F_baratas_dt_train=(2*pred.means2*rec.means2)/(pred.means2+rec.means2)
F_baratas_dt_train
## [1] 0.8775774
tabla_train_arbol
##     pred
## obs    B1   C1
##   B1 7810  924
##   C1 1255 5130
#Media de la F-MEDIDA para Decision Tree
F_dt_train= (F_caras_dt_train+F_baratas_dt_train)/2
F_dt_train
## [1] 0.8512013
accuracy_dt_train = (tabla_train_arbol[1,1]+tabla_train_arbol[2,2]) / (tabla_train_arbol[1,1]+tabla_train_arbol[1,2]+tabla_train_arbol[2,1]+tabla_train_arbol[2,2])
accuracy_dt_train
## [1] 0.8558767

6.3.2 Test

tabla.test.arbol1=table(obs = datos_test_limpio$price_categ1, pred = predict(decisiontree_model,datos_test_limpio[,-15], type ="class"))

tabla.test.arbol1
##     pred
## obs    B1   C1
##   B1 1690  221
##   C1  262 1068
#caras
pred.means2=tabla.test.arbol1[2,2]/(tabla.test.arbol1[2,2]+tabla.test.arbol1[1,2])
rec.means2=tabla.test.arbol1[2,2]/(tabla.test.arbol1[2,2]+tabla.test.arbol1[2,1])
F_medida.means2=(2*pred.means2*rec.means2)/(pred.means2+rec.means2)
F_medida.means2
## [1] 0.8155785
tabla.test.arbol1
##     pred
## obs    B1   C1
##   B1 1690  221
##   C1  262 1068
#baratas
pred.means2=tabla.test.arbol1[1,1]/(tabla.test.arbol1[1,1]+tabla.test.arbol1[2,1])
rec.means2=tabla.test.arbol1[1,1]/(tabla.test.arbol1[1,1]+tabla.test.arbol1[1,2])
F_medida.means2=(2*pred.means2*rec.means2)/(pred.means2+rec.means2)
F_medida.means2
## [1] 0.8749676
tabla.test.arbol1
##     pred
## obs    B1   C1
##   B1 1690  221
##   C1  262 1068
accuracy = (tabla.test.arbol1[1,1]+tabla.test.arbol1[2,2]) / (tabla.test.arbol1[1,1]+tabla.test.arbol1[1,2]+tabla.test.arbol1[2,1]+tabla.test.arbol1[2,2])
accuracy
## [1] 0.8509719

6.3.3 Validacion

tabla.validacion.arbol1=table(obs = datos_validacion_limpio$price_categ1, pred = predict(decisiontree_model,datos_validacion_limpio[,-15], type ="class"))

tabla.validacion.arbol1
##     pred
## obs    B1   C1
##   B1 1700  206
##   C1  263 1074
#caras
pred.means2=tabla.validacion.arbol1[2,2]/(tabla.validacion.arbol1[2,2]+tabla.validacion.arbol1[1,2])
rec.means2=tabla.validacion.arbol1[2,2]/(tabla.validacion.arbol1[2,2]+tabla.validacion.arbol1[2,1])
F_medida.means2=(2*pred.means2*rec.means2)/(pred.means2+rec.means2)
F_medida.means2
## [1] 0.8207872
tabla.validacion.arbol1
##     pred
## obs    B1   C1
##   B1 1700  206
##   C1  263 1074
#baratas
pred.means2=tabla.validacion.arbol1[1,1]/(tabla.validacion.arbol1[1,1]+tabla.validacion.arbol1[2,1])
rec.means2=tabla.validacion.arbol1[1,1]/(tabla.validacion.arbol1[1,1]+tabla.validacion.arbol1[1,2])
F_medida.means2=(2*pred.means2*rec.means2)/(pred.means2+rec.means2)
F_medida.means2
## [1] 0.87878
tabla.validacion.arbol1
##     pred
## obs    B1   C1
##   B1 1700  206
##   C1  263 1074
accuracy = (tabla.validacion.arbol1[1,1]+tabla.validacion.arbol1[2,2]) / (tabla.validacion.arbol1[1,1]+tabla.validacion.arbol1[1,2]+tabla.validacion.arbol1[2,1]+tabla.validacion.arbol1[2,2])
accuracy
## [1] 0.8553808

6.4 Random Forest

6.4.1 Train

randomforest_model=randomForest(price_categ1~., data=datos_train_limpio,parms=list(split="gini") )
#print(randomforest_model)
tabla_randomforest=table(obs = datos_train_limpio$price_categ1, pred = predict(randomforest_model, type = "class") )

tabla_randomforest
##     pred
## obs    B1   C1
##   B1 8001  733
##   C1  726 5659
#caras
pred.means2=tabla_randomforest[2,2]/(tabla_randomforest[2,2]+tabla_randomforest[1,2])
rec.means2=tabla_randomforest[2,2]/(tabla_randomforest[2,2]+tabla_randomforest[2,1])
F_caras_rf_train=(2*pred.means2*rec.means2)/(pred.means2+rec.means2)
F_caras_rf_train
## [1] 0.8858104
tabla_randomforest
##     pred
## obs    B1   C1
##   B1 8001  733
##   C1  726 5659
#baratas
pred.means2=tabla_randomforest[1,1]/(tabla_randomforest[1,1]+tabla_randomforest[2,1])
rec.means2=tabla_randomforest[1,1]/(tabla_randomforest[1,1]+tabla_randomforest[1,2])
F_baratas_rf_train=(2*pred.means2*rec.means2)/(pred.means2+rec.means2)
F_baratas_rf_train
## [1] 0.9164424
tabla_randomforest
##     pred
## obs    B1   C1
##   B1 8001  733
##   C1  726 5659
#Media de la F-MEDIDA para Random Forest
F_rf_train= (F_caras_rf_train+F_baratas_rf_train)/2
F_rf_train
## [1] 0.9011264
accuracy_rf_train = (tabla_randomforest[1,1]+tabla_randomforest[2,2]) / (tabla_randomforest[1,1]+tabla_randomforest[1,2]+tabla_randomforest[2,1]+tabla_randomforest[2,2])
accuracy_rf_train
## [1] 0.9034989

6.4.2 Test

tabla_randomforest_test=table(obs = datos_test_limpio$price_categ1, pred = predict(randomforest_model,datos_test_limpio[,-15], type ="class"))
tabla_randomforest_test
##     pred
## obs    B1   C1
##   B1 1745  166
##   C1  144 1186
#caras
pred.means2=tabla_randomforest_test[2,2]/(tabla_randomforest_test[2,2]+tabla_randomforest_test[1,2])
rec.means2=tabla_randomforest_test[2,2]/(tabla_randomforest_test[2,2]+tabla_randomforest_test[2,1])
F_medida.means2=(2*pred.means2*rec.means2)/(pred.means2+rec.means2)
F_medida.means2
## [1] 0.8844146
tabla_randomforest_test
##     pred
## obs    B1   C1
##   B1 1745  166
##   C1  144 1186
#baratas
pred.means2=tabla_randomforest_test[1,1]/(tabla_randomforest_test[1,1]+tabla_randomforest_test[2,1])
rec.means2=tabla_randomforest_test[1,1]/(tabla_randomforest_test[1,1]+tabla_randomforest_test[1,2])
F_medida.means2=(2*pred.means2*rec.means2)/(pred.means2+rec.means2)
F_medida.means2
## [1] 0.9184211
tabla_randomforest_test
##     pred
## obs    B1   C1
##   B1 1745  166
##   C1  144 1186
accuracy = (tabla_randomforest_test[1,1]+tabla_randomforest_test[2,2]) / (tabla_randomforest_test[1,1]+tabla_randomforest_test[1,2]+tabla_randomforest_test[2,1]+tabla_randomforest_test[2,2])
accuracy
## [1] 0.9043505

6.4.3 Validacion

tabla_randomforest_val=table(obs = datos_validacion_limpio$price_categ1, pred = predict(randomforest_model,datos_validacion_limpio[,-15], type ="class"))

tabla_randomforest_val
##     pred
## obs    B1   C1
##   B1 1756  150
##   C1  148 1189
#caras
pred.means2=tabla_randomforest_val[2,2]/(tabla_randomforest_val[2,2]+tabla_randomforest_val[1,2])
rec.means2=tabla_randomforest_val[2,2]/(tabla_randomforest_val[2,2]+tabla_randomforest_val[2,1])
F_medida.means2=(2*pred.means2*rec.means2)/(pred.means2+rec.means2)
F_medida.means2
## [1] 0.8886398
tabla_randomforest_val
##     pred
## obs    B1   C1
##   B1 1756  150
##   C1  148 1189
#baratas
pred.means2=tabla_randomforest_val[1,1]/(tabla_randomforest_val[1,1]+tabla_randomforest_val[2,1])
rec.means2=tabla_randomforest_val[1,1]/(tabla_randomforest_val[1,1]+tabla_randomforest_val[1,2])
F_medida.means2=(2*pred.means2*rec.means2)/(pred.means2+rec.means2)
F_medida.means2
## [1] 0.9217848
tabla_randomforest_val
##     pred
## obs    B1   C1
##   B1 1756  150
##   C1  148 1189
accuracy = (tabla_randomforest_val[1,1]+tabla_randomforest_val[2,2]) / (tabla_randomforest_val[1,1]+tabla_randomforest_val[1,2]+tabla_randomforest_val[2,1]+tabla_randomforest_val[2,2])
accuracy
## [1] 0.9081098

6.5 SVM

modelo_svm <- e1071::svm(formula = price_categ1 ~., data = datos_train_limpio, kernel = "linear" , probability =TRUE)
summary(modelo_svm)
## 
## Call:
## svm(formula = price_categ1 ~ ., data = datos_train_limpio, kernel = "linear", 
##     probability = TRUE)
## 
## 
## Parameters:
##    SVM-Type:  C-classification 
##  SVM-Kernel:  linear 
##        cost:  1 
## 
## Number of Support Vectors:  4702
## 
##  ( 2352 2350 )
## 
## 
## Number of Classes:  2 
## 
## Levels: 
##  B1 C1
#plot(modelo_svm, datos_train_limpio_dt)
#set.seed(1)
#tune.out = tune(svm,price_categ1 ~ ., data = datos_train_limpio_dt, kernel = "linear", ranges = list(cost = c(0.001,
#0.01, 0.1, 1, 5, 10, 100)))

6.5.1 Train

tabla_svm_train=table(obs = datos_train_limpio$price_categ1, pred = predict(modelo_svm,datos_train_limpio[,-15], type ="class"))

tabla_svm_train
##     pred
## obs    B1   C1
##   B1 7835  899
##   C1 1034 5351
#caras
pred.means2=tabla_svm_train[2,2]/(tabla_svm_train[2,2]+tabla_svm_train[1,2])
rec.means2=tabla_svm_train[2,2]/(tabla_svm_train[2,2]+tabla_svm_train[2,1])
F_caras_svm_train=(2*pred.means2*rec.means2)/(pred.means2+rec.means2)
F_caras_svm_train
## [1] 0.8470123
tabla_svm_train
##     pred
## obs    B1   C1
##   B1 7835  899
##   C1 1034 5351
#baratas
pred_means2=tabla_svm_train[1,1]/(tabla_svm_train[1,1]+tabla_svm_train[2,1])
rec_means2=tabla_svm_train[1,1]/(tabla_svm_train[1,1]+tabla_svm_train[1,2])
F_baratas_svm_train=(2*pred.means2*rec.means2)/(pred.means2+rec.means2)
F_baratas_svm_train
## [1] 0.8470123
tabla_svm_train
##     pred
## obs    B1   C1
##   B1 7835  899
##   C1 1034 5351
#F-Medida para SVM

F_svm_train= (F_caras_svm_train+F_baratas_svm_train)/2
F_svm_train
## [1] 0.8470123
accuracy_svm_train = (tabla_svm_train[1,1]+tabla_svm_train[2,2]) / (tabla_svm_train[1,1]+tabla_svm_train[1,2]+tabla_svm_train[2,1]+tabla_svm_train[2,2])
accuracy_svm_train
## [1] 0.8721476

6.5.2 Test

tabla_svm_test=table(obs = datos_test_limpio$price_categ1, pred = predict(modelo_svm,datos_test_limpio[,-15], type ="class"))
tabla_svm_test
##     pred
## obs    B1   C1
##   B1 1694  217
##   C1  213 1117
#caras
pred.means2=tabla_svm_test[2,2]/(tabla_svm_test[2,2]+tabla_svm_test[1,2])
rec.means2=tabla_svm_test[2,2]/(tabla_svm_test[2,2]+tabla_svm_test[2,1])
F_medida.means2=(2*pred.means2*rec.means2)/(pred.means2+rec.means2)
F_medida.means2
## [1] 0.8385886
tabla_svm_test
##     pred
## obs    B1   C1
##   B1 1694  217
##   C1  213 1117
#baratas
pred.means2=tabla_svm_test[1,1]/(tabla_svm_test[1,1]+tabla_svm_test[2,1])
rec.means2=tabla_svm_test[1,1]/(tabla_svm_test[1,1]+tabla_svm_test[1,2])
F_medida.means2=(2*pred.means2*rec.means2)/(pred.means2+rec.means2)
F_medida.means2
## [1] 0.8873756
tabla_svm_test
##     pred
## obs    B1   C1
##   B1 1694  217
##   C1  213 1117
accuracy = (tabla_svm_test[1,1]+tabla_svm_test[2,2]) / (tabla_svm_test[1,1]+tabla_svm_test[1,2]+tabla_svm_test[2,1]+tabla_svm_test[2,2])
accuracy
## [1] 0.8673249

6.5.3 Validacion

tabla_svm_val=table(obs = datos_validacion_limpio$price_categ1, pred = predict(modelo_svm,datos_validacion_limpio[,-15], type ="class"))

tabla_svm_val
##     pred
## obs    B1   C1
##   B1 1711  195
##   C1  206 1131
#caras
pred.means2=tabla_svm_val[2,2]/(tabla_svm_val[2,2]+tabla_svm_val[1,2])
rec.means2=tabla_svm_val[2,2]/(tabla_svm_val[2,2]+tabla_svm_val[2,1])
F_medida.means2=(2*pred.means2*rec.means2)/(pred.means2+rec.means2)
F_medida.means2
## [1] 0.8494179
#baratas
pred.means2=tabla_svm_val[1,1]/(tabla_svm_val[1,1]+tabla_svm_val[2,1])
rec.means2=tabla_svm_val[1,1]/(tabla_svm_val[1,1]+tabla_svm_val[1,2])
F_medida.means2=(2*pred.means2*rec.means2)/(pred.means2+rec.means2)
F_medida.means2
## [1] 0.8951086
tabla_svm_val
##     pred
## obs    B1   C1
##   B1 1711  195
##   C1  206 1131
accuracy = (tabla_svm_val[1,1]+tabla_svm_val[2,2]) / (tabla_svm_val[1,1]+tabla_svm_val[1,2]+tabla_svm_val[2,1]+tabla_svm_val[2,2])
accuracy
## [1] 0.8763491

7 GAM

7.0.1 Train

datos_train_gam <- datos_train_limpio
datos_train_gam$price <- datos_train$price
 
model_gam <- gam(price ~ s(bedrooms, bs = "cr") + bathrooms_group + s(log_sqft_living, bs = "cr") + 
                   s(log_lot, bs = "cr") + sqft_basement_cat + waterfront + view + s(lat, bs = "cr") + 
                   s(long, bs = "cr") + zona + yr_renovated_catg, data = datos_train_gam)

summary(model_gam)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## price ~ s(bedrooms, bs = "cr") + bathrooms_group + s(log_sqft_living, 
##     bs = "cr") + s(log_lot, bs = "cr") + sqft_basement_cat + 
##     waterfront + view + s(lat, bs = "cr") + s(long, bs = "cr") + 
##     zona + yr_renovated_catg
## 
## Parametric coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          460094       3772 121.978  < 2e-16 ***
## bathrooms_group1      35433       4091   8.662  < 2e-16 ***
## sqft_basement_cat1   -33249       3315 -10.029  < 2e-16 ***
## waterfront1          523178      20514  25.503  < 2e-16 ***
## view1                 90905      11827   7.686 1.61e-14 ***
## view2                 82058       7077  11.595  < 2e-16 ***
## view3                155689       9712  16.031  < 2e-16 ***
## view4                287662      14672  19.606  < 2e-16 ***
## zona2                126977       5000  25.398  < 2e-16 ***
## yr_renovated_catg1    58719       7070   8.306  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                      edf Ref.df      F p-value    
## s(bedrooms)        8.739  8.957   19.8  <2e-16 ***
## s(log_sqft_living) 8.911  8.997 1509.6  <2e-16 ***
## s(log_lot)         8.339  8.802   37.3  <2e-16 ***
## s(lat)             8.872  8.995  280.6  <2e-16 ***
## s(long)            8.935  8.999  171.2  <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.781   Deviance explained = 78.2%
## GCV = 3.0207e+10  Scale est. = 3.0099e+10  n = 15119
plot(model_gam, ylab = "price")

par(mfrow=c(3,3))
visreg(model_gam)

7.0.2 Test

datos_test_limpio$price <- datos_test$price
datos_test_limpio$price_predict = predict(model_gam, datos_test_limpio)


datos_test_limpio$precio_diff <- abs(datos_test_limpio$price - datos_test_limpio$price_predict)

dif_media_test <- mean(datos_test_limpio$precio_diff)
dif_media_test
## [1] 108157.9

7.0.3 Validación

datos_validacion_limpio$price <- datos_validacion$price
datos_validacion_limpio$price_predict = predict(model_gam, datos_validacion_limpio)


datos_validacion_limpio$precio_diff <- abs(datos_validacion_limpio$price - datos_validacion_limpio$price_predict)

dif_media_test <- mean(datos_validacion_limpio$precio_diff)
dif_media_test
## [1] 105924.8

8 Ajuste de Hiperparámetros.

#knn.train2=knn.cv(scale(datos_train_numeric),cl=as.factor(datos_train_limpio$pri#ce_categ1),k=2)

 # x <- iris[,-5]
  #y <- iris[,5]
  
obj2 <- tune.knn(x=scale(datos_train_numeric),
                   y=as.factor(datos_train_limpio$price_categ1), k = 1:10,
                   tunecontrol = tune.control(sampling = "boot"))
 #summary(obj2)
  plot(obj2)

8.1 Ajuste KNN

#n^(4/(d+4))
#n es el número de datos de entrenamiento
#d es la dimensión de los datos
n=dim(datos_train_numeric)[1]
d=ncol(datos_train_numeric)
k_optimo = as.integer(n^(4/(d+4)))

k_optimo
## [1] 72
rango=1:(k_optimo+10)
modelos=c()
f1_modelos=c()
for (i in rango){
  model_knn=knn.cv(scale(datos_train_numeric),cl=as.factor(datos_train_limpio$price_categ1),k=i)
  tabla=table(datos_train_limpio$price_categ1,model_knn)
  # f1 casas caras
  pred_means_caras=tabla[2,2]/(tabla[2,2]+tabla[1,2])
  rec_means_caras=tabla[2,2]/(tabla[2,2]+tabla[2,1])
  f1_caras=(2*pred_means_caras*rec_means_caras)/(pred_means_caras+rec_means_caras)
  # f1 casas baratas
  pred_means_baratas=tabla[1,1]/(tabla[1,1]+tabla[2,1])
  rec_means_baratas=tabla[1,1]/(tabla[1,1]+tabla[1,2])
  f1_baratas=(2*pred_means_baratas*rec_means_baratas)/(pred_means_baratas+rec_means_baratas)
  f1_total = (f1_baratas + f1_caras)/2
  f1_modelos=c(f1_modelos,f1_total)
}


f1_modelos
##  [1] 0.8575267 0.8503791 0.8735162 0.8683478 0.8753198 0.8755879 0.8766022
##  [8] 0.8751717 0.8786105 0.8752529 0.8780367 0.8790804 0.8792161 0.8771628
## [15] 0.8787699 0.8775224 0.8797486 0.8789995 0.8782733 0.8796625 0.8784356
## [22] 0.8782627 0.8793153 0.8790047 0.8790228 0.8788375 0.8785713 0.8780670
## [29] 0.8787176 0.8775505 0.8775171 0.8767733 0.8767294 0.8764420 0.8769355
## [36] 0.8766922 0.8764966 0.8755295 0.8766043 0.8753725 0.8760616 0.8751011
## [43] 0.8745799 0.8743301 0.8743193 0.8731636 0.8740480 0.8743576 0.8743301
## [50] 0.8736573 0.8742159 0.8732341 0.8733860 0.8729248 0.8728759 0.8729464
## [57] 0.8729518 0.8716501 0.8719591 0.8709986 0.8711076 0.8710807 0.8706303
## [64] 0.8711186 0.8708801 0.8710535 0.8701859 0.8700173 0.8692036 0.8699356
## [71] 0.8698332 0.8701804 0.8694319 0.8698662 0.8689657 0.8694485 0.8691829
## [78] 0.8696820 0.8687928 0.8678551 0.8685051 0.8676693
rango
##  [1]  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
## [26] 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50
## [51] 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75
## [76] 76 77 78 79 80 81 82
pos_max <- which.max(f1_modelos)
rango[pos_max]
## [1] 17

9 Evaluación y comparación de modelos.

Una vez que se han entrenado y optimizado distintos modelos, se tiene que identificar cuál de ellos consigue mejores resultados para el problema en cuestión, en este caso, predecir si una casa es barata o cara. Con los datos disponibles, existen dos formas de comparar los modelos. Si bien las dos no tienen por qué dar los mismos resultados, son complementarias a la hora de tomar una decisión final.

models_cross = data.frame(
"modelo"= c('GLM','knn','Decision_Tree','Random_Forest','SVM'),
 "F_Medida_train" = c(F_reglog_train,F_knn_train,F_dt_train,F_rf_train,F_svm_train))
plot(x = models_cross$modelo, y= models_cross$F_Medida_train,fill=models_cross$modelo)

ggplot(data=models_cross, aes(x=modelo, y=F_Medida_train, fill=modelo)) + 
    geom_bar(stat="identity", position="dodge")

10 CURVA ROC.

# REGRESIÓN LOGÍSTICA
predictions_glm <- predict(model_glm, newdata = datos_train_rl, type = "response")
pred_glm <- prediction(predictions_glm, datos_train_rl$price_categ1)
perf_glm <- performance(pred_glm,"tpr","fpr")

# ARBOL DE DECISION
predictions_tree <- predict(decisiontree_model, newdata = datos_decision_tree, type = "prob")
pred_tree = prediction(predictions_tree[,2], datos_decision_tree$price_categ1)
pref_tree = performance(pred_tree, "tpr", "fpr")

# RANDOM FOREST
predictions_rf <- predict(randomforest_model, new_data=datos_train_limpio, type = "prob")
pred_rf <- prediction(predictions_rf[,2],datos_train_limpio$price_categ1)
perf_rf <- performance(pred_rf,"tpr","fpr")

# SVM
predictions_svm <- predict(modelo_svm, newdata=datos_train_limpio, probability = TRUE)
prob_svm<-attr(predictions_svm,"probabilities")
pred_svm <- prediction(prob_svm[,2], datos_train_limpio$price_categ)
perf_svm <- performance(pred_svm,"tpr","fpr")


plot(perf_glm,col="blue")
plot(pref_tree, col="red",add = TRUE)
plot(perf_rf,col="green", add = TRUE)
plot(perf_svm, col="yellow",add = TRUE)
legend(x="right" ,legend=c("GLM","DT","RF","SVM"),
fill=c("blue","red","green","yellow"),cex=0.8)

11 Ideas

cambiar el modelo para que aumente 1000 más la variable respuesta

11.1 para estructurarlo

Primero transformaciones de train, test validacion Luego cada modelo por separado con su train y test y validacion

11.2 objetivo

11.3 concluisones

  • nos quedamos con la acuracy ya que es la única medida que tiene en cuenta los aciertos de las casas caras y de las baratas (la F_medida no tiene en cuenta los TN y la sensibilidad y especificidad solo tienen en cuenta que acierte las caras o las baratas pero no que acierte en ambas)